
#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include "Header.h"
#include <stdio.h>

char *cc;

__device__ void d_trans_refer_to_digits(const char* const refer_seq, 
										const int* const SA, 
										const char* const trans,
										const int thx,
										const int index,
										const int n,
										const int c,
										unsigned int* s_ref)
{
	for(int i = 0; i < 4; i++) //16bp characters -> 16*2(bits) = 32 bits = (integer type)4 bytes -> 128bp / 16bp = 8
	{
 		unsigned int result = 0;
		int t = 30;
		for(int j = 0; j < 16; j++)
		{
			if(c * 64 + SA[index] + i*16 + j < n) 
				result += (((unsigned int)trans[ refer_seq[SA[index] + c * 64 + i*16 + j] - 'A']) << t); //transform 16bp to an integer
			t -= 2;
		}
		s_ref[thx * 4 + i] = result;
	}

}

__device__ int d_query_refer_cmp(int thx, int c, unsigned int s_ref[], unsigned int s_query[])
{
	for(int i = 0; i < 4; i++)
	{
		if (s_query[thx % 128 + 4 * 128 * c + 128 * i] != s_ref[thx * 4 + i])
			return s_query[thx % 128 + 4 * 128 * c + 128 * i] > s_ref[thx * 4 + i] ? 1 : -1;
	}
	return 0;
	
}

__global__ void gpu_genome_matching(const int* const SA, 
									const int n, 
									unsigned int* query_seq, 
									const char* const refer_seq, 
									int* boundaries,
									char* trans)
{
	
	int thx =  blockDim.x * blockIdx.x + threadIdx.x;
	int boundary = 0, t, i, j, k,temp;

	__shared__ unsigned int s_query[128*8], s_ref[256*4];


	for(i = 0; i < 8; i++)
	{
		s_query[thx % 128 + 128 * i] = query_seq[thx % 128 * 8 + i];
	}
	//0010 0111 1111 1101 0011 1101 1110 0001
    __syncthreads();

		if (thx < 128) // find LB
		{

			for(t = 0,i = 0 ; i < 2; i++)
			{		
				//d_trans_refer_to_digits(refer_seq, SA, trans, thx, 0, n,i, s_ref);
				for( k = 0; k < 4; k++) //16bp characters -> 16*2(bits) = 32 bits = (integer type)4 bytes -> 128bp / 16bp = 8
				{ 
					unsigned int result = 0;
					temp = 30;
					for( j = 0; j < 16; j++)
					{
						if(SA[0] + i * 64 + k*16 + j < n) 
							result += (((unsigned int)trans[ refer_seq[SA[0] + i * 64 + k*16 + j] - 'A']) << temp); //transform 16bp to an integer
						temp -= 2;
					}
					s_ref[thx * 4 + k] = result;
				}
				 
				for(j = 0; j < 4; j++)
				{
					if (s_query[thx + 4 * 128 * i + 128 * j] != s_ref[thx * 4 + j])
					{
						t = s_query[thx + 4 * 128 * i + 128 * j] > s_ref[thx * 4 + j] ? 1 : -1;
						break;
					}
				
				}

				//t += d_query_refer_cmp(thx,i,s_ref,s_query);

				if(t != 0) 
				{
					if(t < 0)
						boundary = t;
					break;
				}
			}
			if (t == 1)
			{

				for(t = 0,i = 0 ; i < 2; i++)
				{		
					//d_trans_refer_to_digits(refer_seq, SA, trans, thx, n-1, n,i, s_ref);
			
					for( k = 0; k < 4; k++) //16bp characters -> 16*2(bits) = 32 bits = (integer type)4 bytes -> 128bp / 16bp = 8
					{
						unsigned int result = 0;
						temp = 30;
						for( j = 0; j < 16; j++)
						{
							if(SA[n-1] + i * 64 + k*16 + j < n) 
								result += (((unsigned int)trans[refer_seq[SA[n-1] + i * 64 + k * 16 + j] - 'A']) << temp); //transform 16bp to an integer
							temp -= 2;
						}
						s_ref[thx * 4 + k] = result;
					}
					//__syncthreads();
					//t += d_query_refer_cmp(thx,i,s_ref,s_query);
					for(j = 0; j < 4; j++)
					{

						if (s_query[thx + 4 * 128 * i + 128 * j] != s_ref[thx * 4 + j])
						{
							t = s_query[thx + 4 * 128 * i + 128 * j] > s_ref[thx * 4 + j] ? 1 : -1;
							break;
						}
					}


					if(t != 0) 
					{
						if(t > 0)
							boundary = -t;
						break;
					}
				}
				if (t != 1)
				{
					int m, l = 0, r = n;
					while(r > l + 1)
					{
						m = (l + r)/2;

						for(t = 0,i = 0 ; i < 2; i++)
						{		
							//d_trans_refer_to_digits(refer_seq, SA, trans, thx, m, n,i, s_ref);
						
							for( k = 0; k < 4; k++) //16bp characters -> 16*2(bits) = 32 bits = (integer type)4 bytes -> 128bp / 16bp = 8
							{
								unsigned int result = 0;
								temp = 30;
								for( j = 0; j < 16; j++)
								{
									if(SA[m] + i * 64 + k*16 + j < n) 
										result += (((unsigned int)trans[refer_seq[SA[m] + i * 64 + k * 16 + j] - 'A']) << temp); //transform 16bp to an integer
									temp -= 2;
								}
								s_ref[thx * 4 + k] = result;
							}

							//t += d_query_refer_cmp(thx,i,s_ref,s_query);
							for(j = 0; j < 4; j++)
							{
								if (s_query[thx + 4 * 128 * i + 128 * j] != s_ref[thx * 4 + j])
								{
									t = s_query[thx + 4 * 128 * i + 128 * j] > s_ref[thx * 4 + j] ? 1 : -1;
									break;
								}
							}

							if(t != 0) 
								break;
						}
						if(t <= 0)
							r = m;
						else 
							l = m;
					}
					boundary = r;
				}

			}

		}else  // find RB
		{
			for(t = 0,i = 0 ; i < 2; i++)
			{
				//d_trans_refer_to_digits(refer_seq, SA, trans, thx, 0, n,i, s_ref);

				for( k = 0; k < 4; k++) //16bp characters -> 16*2(bits) = 32 bits = (integer type)4 bytes -> 128bp / 16bp = 8
				{
					unsigned int result = 0;
					temp = 30;
					for( j = 0; j < 16; j++)
					{
						if(SA[0] + i * 64 + k*16 + j < n)  
							result += (((unsigned int)trans[refer_seq[SA[0] + i * 64 + k * 16 + j] - 'A']) << temp); //transform 16bp to an integer
						temp -= 2;
					}
					s_ref[thx * 4 + k] = result;
				}

				//t += d_query_refer_cmp(thx,i,s_ref,s_query);
				for(j = 0; j < 4; j++)
				{
					if (s_query[thx % 128 + 4 * 128 * i + 128 * j] != s_ref[thx * 4 + j])
					{
						t = s_query[thx % 128 + 4 * 128 * i + 128 * j] > s_ref[thx * 4 + j] ? 1 : -1;
						break;
					}
				}
				
				
				if(t != 0) 
				{
					if(t < 0)
						boundary = t;
					break;
				}
			}
			if (t >= 0)
			{

				for(t = 0,i = 0 ; i < 2; i++)
				{		
					//d_trans_refer_to_digits(refer_seq, SA, trans, thx, n-1, n,i,s_ref);

					for( k = 0; k < 4; k++) //16bp characters -> 16*2(bits) = 32 bits = (integer type)4 bytes -> 128bp / 16bp = 8
					{
						unsigned int result = 0;
						temp = 30;
						for( j = 0; j < 16; j++)
						{
							if(SA[n - 1] + i * 64 + k*16 + j < n) 
								result += (((unsigned int)trans[refer_seq[SA[n - 1] + i * 64 + k * 16 + j] - 'A']) << temp); //transform 16bp to an integer
							temp -= 2;
						}
						s_ref[thx * 4 + k] = result;
					}

					//t += d_query_refer_cmp(thx,i,s_ref,s_query);
					for(j = 0; j < 4; j++)
					{
						if (s_query[thx % 128 + 4 * 128 * i + 128 * j] != s_ref[thx * 4 + j])
						{
							t = s_query[thx % 128 + 4 * 128 * i + 128 * j] > s_ref[thx * 4 + j] ? 1 : -1;
							break;
						}
					}

					if(t != 0) 
					{
						if(t > 0)
							boundary = -t;
						break;
					}
				}
				if (t == 0)
					boundary = n - 1;
				else if (t == -1)
				{
					int m, l = 0, r = n;
					while(r > l + 1)
					{
						m = (l + r)/2;

						for(t = 0,i = 0 ; i < 2; i++)
						{		
							
							
							//d_trans_refer_to_digits(refer_seq, SA, trans, thx, m, n,i,s_ref);
							for( k = 0; k < 4; k++) //16bp characters -> 16*2(bits) = 32 bits = (integer type)4 bytes -> 128bp / 16bp = 8
							{
								unsigned int result = 0;
								temp = 30;
								for( j = 0; j < 16; j++)
								{
									if(SA[m] + i * 64 + k*16 + j < n)  
										result += (((unsigned int)trans[refer_seq[SA[m] + i * 64 + k * 16 + j] - 'A']) << temp); //transform 16bp to an integer
									temp -= 2;
								}
								s_ref[thx * 4 + k] = result;
							}

							//t += d_query_refer_cmp(thx,i,s_ref,s_query);
							for(j = 0; j < 4; j++)
							{
								if (s_query[thx % 128 + 4 * 128 * i + 128 * j] != s_ref[thx * 4 + j])
								{
									t = s_query[thx % 128 + 4 * 128 * i + 128 * j] > s_ref[thx * 4 + j] ? 1 : -1;
									break;
								}
							}
							if(t != 0) 
								break;
						}
						if(t >= 0)
							l = m;
						else 
							r = m;
					}
					boundary = l;
				}	
			}
		}
		//printf("thx:%d B: %d", thx, boundary);
		boundaries[thx] = boundary;
}

int main(int argc, char* argv[])
{
	//freopen("genome","r",stdin);
	//freopen("output.txt","w",stdout);
	clock_t start, end;
	double runTime;
	char *filename = "genome";
	int num = 10000;//modify the length of genome data 
	char *genome;
	int n = 0,i = 0;
	int *inp;
	int *SA;
	genome = (char *) malloc((num+1)*sizeof(char)); 

	read_genome(filename, genome, num);

	start = clock();
	cc = genome;
	n = strlen(cc);

	inp = (int *)malloc( (n+3)*sizeof(int) );
	SA  = (int *)malloc( (n+3)*sizeof(int) );

	for(i=0;i<n;i++)
		inp[i] = to_i(cc[i]);
	inp[i]=0;inp[i+1]=0;inp[i+2]=0;

	for(i=0;i<n+3;i++)
		SA[i] = 0;

	suffixArray(inp,SA,n,MAX_ALPHA);

	//  print all suffixes of genome
	//	for(i = 0 ; i < num ; i++)
	//		print_suffix(cc, SA[i]);

	end = clock();
	runTime = (end - start) / (double) CLOCKS_PER_SEC ;

	printf("CPU linear construct Suffix Array\nNUM: %d \t Time: %f Sec\n", num, runTime);

	//on CPU, define LB and RB 
	/*
	int LB = CPU_left_boundary(SA, num, "AT", cc);
	int RB;
	if(LB != -1)
	{
		RB = CPU_right_boundary(SA, num, "ATC", cc);

		//  print all indexed suffixes of genome
		//	for(i = LB ; i <= RB ; i++)
		//		print_suffix(cc, SA[i]);
	}
	*/
	unsigned int h_query[128][8];
	for(i = 0; i < 128; i++)
	{
		int t = i * 10;
		if(i == 3)
			t = 174;
		trans_querySeq_to_digits(&cc[t], h_query[i]); //random generate #256 reference sequences(each has 1024bp) 
	}
	
	int h_left_right_boundary[128*2]; //a thread block processes 128 query sequences and results need 128*2 (RB and LB)
	int* d_left_right_boundary;
	unsigned int* d_query;
	char* d_genome;
	int*  d_SA;
	dim3 gridDim(1,1,1);
	dim3 blockDim(256,1,1);
	char trans[MAX_ALPHA];
	memset(trans, 0, MAX_ALPHA);
	char* d_trans;
	trans[0] = 0;  // 00 'A'- 'A'
	trans[2] = 1;  // 01 'C'- 'A'
	trans[6] = 2;  // 10 'G'- 'A'
	trans[19] = 3;  // 11 'T'- 'A'

	cudaMalloc((void**)&d_left_right_boundary, sizeof(int)*128*2);
	cudaMalloc((void**)&d_query, sizeof(unsigned int)*128*8);
	cudaMalloc((void**)&d_genome, sizeof(char)*n);
	cudaMalloc((void**)&d_SA, sizeof(int)*n);
	cudaMalloc((void**)&d_trans, sizeof(char)*MAX_ALPHA);

	cudaMemcpy(d_trans, trans, sizeof(char)*MAX_ALPHA, cudaMemcpyHostToDevice);
	cudaMemcpy(d_query, h_query, sizeof(unsigned int)*128*8, cudaMemcpyHostToDevice);
	cudaMemcpy(d_genome, genome, sizeof(char)*n, cudaMemcpyHostToDevice);
	cudaMemcpy(d_SA, SA, sizeof(int)*n, cudaMemcpyHostToDevice);

	gpu_genome_matching<<<gridDim, blockDim>>>(d_SA, n, d_query, d_genome, d_left_right_boundary, d_trans);
	
	cudaMemcpy(h_left_right_boundary, d_left_right_boundary, sizeof(int)*128*2, cudaMemcpyDeviceToHost);

	for (i = 0; i < 128; i++)
	{
		
		printf("%d %d %d\n", i, h_left_right_boundary[i], h_left_right_boundary[i + 128]);

		for(int j = h_left_right_boundary[i]; j <= h_left_right_boundary[i + 128]; j++)
		{
			int k = 0;
			printf("%d: ",j);
			for(k = SA[j]; k < SA[j] + 128; k++)
				printf("%c",cc[k]);
			printf("\n");
		}
	}
	
	
	
	free(genome);
	return 0;
}



/*
int CPU_left_boundary(int* SA, int n, char* query_seq, char* refer_seq)
{
	int m, temp;
	if (strncmp(query_seq,&refer_seq[SA[0]],strlen(query_seq)) == 0)
		return 0;

	if(strncmp(query_seq,&refer_seq[SA[0]],strlen(query_seq)) < 0)
		return -1;
	else if(strncmp(query_seq,&refer_seq[SA[n-1]],strlen(query_seq)) > 0)
		return -1;
	else
	{
		int l = 0;
		int r = n;
		while(r > l+1)
		{
			m = (l + r)/2;
			temp = strncmp(query_seq,&refer_seq[SA[m]],strlen(query_seq));

			if(temp <= 0)
				r = m;
			else 
				l = m;
		}
		return r;
	}
}

int CPU_right_boundary(int* SA, int n, char* query_seq, char* refer_seq)
{
	if(strncmp(query_seq,&refer_seq[SA[n-1]],strlen(query_seq)) == 0)
		return n-1;
	int m, temp;

	int l = 0;
	int r = n;
	while(r > l + 1 )
	{
		m = (l + r)/2;
		temp = strncmp(query_seq,&refer_seq[SA[m]],strlen(query_seq));

		if(temp >= 0)
			l = m;
		else 
			r = m;
	}
	return l;

}*/